Import data

# Import data
df <- readr::read_rds('./data/Juliane_raw_FEST.rds')

Quick look

# Quick look
dim(df)
## [1] 2735   11
head(df)
## # A tibble: 6 x 11
## # Groups:   id [1]
##      id block_number trial_no trial_type   perc_trial_type CS      US   
##   <dbl> <chr>           <int> <chr>        <chr>           <chr>   <chr>
## 1  8.00 baseline            1 CS_only      CS_only         CSplus  <NA> 
## 2  8.00 baseline            2 CSminus_only CSminus_only    CSminus <NA> 
## 3  8.00 baseline            3 CS_only      CS_only         CSplus  <NA> 
## 4  8.00 baseline            4 CS_only      CS_only         CSplus  <NA> 
## 5  8.00 baseline            5 CSminus_only CSminus_only    CSminus <NA> 
## 6  8.00 baseline            6 CSminus_only CSminus_only    CSminus <NA> 
## # ... with 4 more variables: rating <int>, rating_time <int>,
## #   button_correct <chr>, button_time <int>
tail(df)
## # A tibble: 6 x 11
## # Groups:   id [1]
##      id block_number trial_no trial_type    perc_trial_type CS      US   
##   <dbl> <chr>           <int> <chr>         <chr>           <chr>   <chr>
## 1  30.0 test3             137 CSminus_USlow CSminus_USlow   CSminus low  
## 2  30.0 test3             138 CSminus_only  CSminus_only    CSminus <NA> 
## 3  30.0 test3             139 CSplus_UShigh CSplus_UShigh   CSplus  high 
## 4  30.0 test3             140 CSminus_only  CSminus_only    CSminus <NA> 
## 5  30.0 test3             141 CSminus_USlow CSminus_USlow   CSminus low  
## 6  30.0 test3             142 CSplus_UShigh CSplus_UShigh   CSplus  high 
## # ... with 4 more variables: rating <int>, rating_time <int>,
## #   button_correct <chr>, button_time <int>
glimpse(df)
## Observations: 2,735
## Variables: 11
## $ id              <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8...
## $ block_number    <chr> "baseline", "baseline", "baseline", "baseline"...
## $ trial_no        <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,...
## $ trial_type      <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ perc_trial_type <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ CS              <chr> "CSplus", "CSminus", "CSplus", "CSplus", "CSmi...
## $ US              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "low",...
## $ rating          <int> 0, -7, -7, -24, -32, -39, -17, -33, -40, -33, ...
## $ rating_time     <int> 4824, 4342, 5833, 4173, 2689, 3275, 3460, 3201...
## $ button_correct  <chr> "correct", "correct", "correct", "correct", "c...
## $ button_time     <int> 8407, 5485, 7576, 5572, 3792, 4610, 4875, 4152...

Plot ratings by trial type

Ratings on the FEST are plotted for each participant and then for the whole sample. First we plot by trial type delivered, then by trial type perceived.

# Make plots
ggplot(data = df) +
  aes(x = trial_no,
      y = rating,
      colour = trial_type) +
  geom_point() +
  scale_colour_brewer(palette = "Set1") +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 10.5, linetype = 3) +
  geom_vline(xintercept = 34.5, linetype = 3) +
  facet_wrap(~ id) +
  labs(title = "Ratings over course of procedure",
       subtitle = "Faceted by ID")

df_test <- df %>% filter(trial_type == "CSplus_UStest" | trial_type == "CSminus_UStest")

ggplot(data = df_test) +
  aes(x = trial_no,
      y = rating,
      colour = trial_type) +
  geom_point() +
  scale_colour_brewer(palette = "Set1") +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 70.5, linetype = 3) +
  geom_vline(xintercept = 34.5, linetype = 3) +
  geom_vline(xintercept = 105.5, linetype = 3) +
  facet_wrap(~ id) +
  labs(title = "Ratings over course of test phases, UStest trials only",
       subtitle = "Faceted by ID")

ggplot(data = df_test) +
  aes(x = as.factor(trial_type),
      y = rating,
      colour = trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  facet_wrap(~ id) +
  labs(title = "Ratings of test stimuli according to actual trial type",
       subtitle = "Faceted by ID")

ggplot(data = df_test) +
  aes(x = as.factor(perc_trial_type),
      y = rating,
      colour = perc_trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  facet_wrap(~ id) +
  labs(title = "Ratings of test stimuli according to perceived trial type",
       subtitle = "Faceted by ID")

df_test %>% filter(!is.na(trial_type)) %>%
ggplot(data = .) +
  aes(x = as.factor(trial_type),
      y = rating,
      colour = trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Ratings of test stimuli",
       subtitle = "According to actual trial type")

df_test %>% filter(!is.na(perc_trial_type)) %>%
  ggplot(data = .) +
  aes(x = as.factor(perc_trial_type),
      y = rating,
      colour = perc_trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Ratings of test stimuli",
       subtitle = "According to perceived trial type")

Prediction interval plots

Here, we plot the prediction interval for the ‘control’ stimulua type (e.g. CSminus/UStest) and depict the trials of the comparator trial type (e.g. CSplus/UStest) that fall within and outside the predictio interval. We do this for the test trials, as delivered and as perceived, and then for the acquisition trials.

# Prediction intervals for test trials as DELIVERED
df_prediction <- df_test %>%
  # Filter for CSminus_UStest
  filter(trial_type == 'CSminus_UStest' & id != 23) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df_test) %>%
  # Is CSplus_UStest within the prediction interval
  mutate(CC_effect = ifelse(rating > upper_pi,
                               yes = 'yes',
                               no = 'no')) %>%
  # Filter for CSplus_UStest
  filter(trial_type == 'CSplus_UStest') 

ggplot(data = df_prediction) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = CC_effect),
             height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for CSminus_UStest trial ratings',
       subtitle = 'Points are individual CSplus_UStest trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

# Prediction intervals for test trials as PERCEIVED

id_19 <- df_test %>% filter(id == 19)

plot(id_19$rating, id_19$trial_no)

df_prediction_perc <- df_test %>%
  filter(!is.na(perc_trial_type)) %>%
  # Filter for CSminus_UStest
  filter(perc_trial_type == 'CSminus_UStest' & id != 23 & id != 19) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df_test) %>%
  # Is CSplus_UStest within the prediction interval
  mutate(CC_effect = ifelse(rating > upper_pi,
                            yes = 'yes',
                            no = 'no')) %>%
  # Filter for CSplus_UStest
  filter(perc_trial_type == 'CSplus_UStest') 

# Save file df_prediction_perc

write_rds(x = df_prediction_perc,
          path = './data/Juliane_raw_as_perc.rds')

ggplot(data = df_prediction_perc) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = CC_effect),
              height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for ratings of trials perceived as CSminus_UStest ',
       subtitle = 'Points are individual CSplus_UStest (percieved) trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

# Prediction intervals for all reinforced (CS+/USH or CS-/USL) trials

df_prediction_reinf <- df %>%
  # Filter for CSminus_USlow
  filter(trial_type == 'CSminus_USlow' & id != 23) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df) %>%
  # Is CSplus_UShigh within the prediction interval
  mutate(diff_intensities = ifelse(rating > upper_pi,
                            yes = 'yes',
                            no = 'no')) %>%
  # Filter for CSplus_UShigh
  filter(trial_type == 'CSplus_UShigh') 

ggplot(data = df_prediction_reinf) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = diff_intensities),
              height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for CSminus_USlow trial ratings',
       subtitle = 'Points are individual CSplus_UShigh trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

Conclusion about calibration

Calibration was poor. However, participants 11, 15 and 26 have a prediction interval for their CS-/USlow trials that does not cross zero, and they also show no CC effect according to the previous plot.

Session information

sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2       RColorBrewer_1.1-2 magrittr_1.5      
##  [4] dplyr_0.7.4        purrr_0.2.4        tidyr_0.7.2       
##  [7] tibble_1.4.2       ggplot2_2.2.1.9000 tidyverse_1.1.1   
## [10] readr_1.1.1       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.15      cellranger_1.1.0  pillar_1.1.0     
##  [4] compiler_3.4.3    plyr_1.8.4        bindr_0.1        
##  [7] forcats_0.2.0     tools_3.4.3       digest_0.6.15    
## [10] lubridate_1.6.0   jsonlite_1.5      evaluate_0.10.1  
## [13] nlme_3.1-131      gtable_0.2.0      lattice_0.20-35  
## [16] pkgconfig_2.0.1   rlang_0.1.6.9003  psych_1.7.8      
## [19] cli_1.0.0         yaml_2.1.14       parallel_3.4.3   
## [22] haven_1.1.0       withr_2.1.1.9000  xml2_1.1.1       
## [25] httr_1.3.1        stringr_1.2.0     knitr_1.17       
## [28] hms_0.3           rprojroot_1.2     grid_3.4.3       
## [31] glue_1.1.1        R6_2.2.2          readxl_1.0.0     
## [34] foreign_0.8-69    rmarkdown_1.6     modelr_0.1.1     
## [37] reshape2_1.4.3    backports_1.1.1   scales_0.5.0.9000
## [40] htmltools_0.3.6   rvest_0.3.2       assertthat_0.2.0 
## [43] mnormt_1.5-5      colorspace_1.3-2  labeling_0.3     
## [46] utf8_1.1.3        stringi_1.1.5     lazyeval_0.2.1   
## [49] munsell_0.4.3     broom_0.4.2       crayon_1.3.4